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Unsteady  Laminar  Flow  in  a Circular  Tube: 

A test  of  the  HERCOL  (Hermitian  collocation)  Computer  Code 

James  F.  Welch,  James  A.  Hurley,  Michael  P.  Glover, 
Ryan  D.  Nassimbene,  and  Marilyn  R.  Yetzbacher 


HERCOL,  a computer  code  for  the  integration  of  second— order  differential  equations  in  one 
space  dimension  by  Hermitian  collocation,  was  used  to  calculate  the  unsteady  velocity  profiles  for  laminar 
flow  in  a circular  tube.  The  code  was  tested  for  stability  and  accuracy  on  this  problem  for  which  an 
analytical  solution  exists  prior  to  application  to  a like  problem  in  which  the  initial  and  boundary 
conditions  preclude  the  existence  of  analytical  solutions. 

The  test  problem  is  one  in  which  a pressure  gradient  is  imposed  on  a fluid  initially  at  rest  in  a 
circular  tube;  the  fluid  accelerates,  and,  at  steady  state,  has  a parabolic  velocity  profile.  A second 
example  was  constructed  from  the  first;  a pressure  gradient  equal  but  opposite  in  sign  is  imposed  on  the 
fluid  with  a fully  developed  parabolic  velocity  profile.  At  steady  state,  the  velocity  is  again  parabolic  but 
in  the  opposite  direction  to  that  at  the  initial  conditions. 

Excellent  agreement  with  the  analytical  solution  was  obtained  in  the  first  problem;  in  the  second, 
the  behaviour  was  as  expected.  This  example  is  suitable  for  first— time  users  of  the  code. 

Key  words:  numerical  integration;. partial  differential  equation;  unsteady— state  laminar  flow; 
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Unsteady— State  Laminar  Velocity  Profiles 


INTRODUCTION 

Velocity  profiles  for  two  examples  of  unsteady  laminar  flow  of  a fluid  of  constant  density,  /?,  and 
constant  viscosity,  //,  in  a horizontal  tube  of  length  L and  radius  R (Figure  1)  were  calculated  with  the 
aid  of  HERCOL  [1],  a computer  code  for  the  integration  of  second— order  differential  operators  in  one 
space  dimension.  For  the  first  example,  the  analytical  solution  [2,3]  was  compared  with  the  results 
obtained  from  the  code  to  assess  the  performance  of  HERCOL  and  to  estimate  its  efficacy  for  the 
solution  of  similar  problems  for  which  no  analytical  solution  exists.  A simple  modification  of  the  first 
example  to  form  a second  was  used  to  emphasize  the  speed  and  ease  with  which  problems,  with  no 
readily  available  analytical  solution  may  be  solved  with  HERCOL. 

In  addition  to  serving  as  a test  for  HERCOL,  the  solution  may  be  useful  for  estimating  the  time 
required  for  the  flow  in  a sampling  tube  operated  in  a periodic  fashion  to  achieve  steady  state. 

UNSTEADY  LAMINAR  FLOW  IN  A CIRCULAR  TUBE 

A pressure  gradient,  — of  constant  magnitude  (p^  — imposed  on  a fluid  contained 

in  a tube  shown  in  Figure  1.  The  fluid  initially  at  rest  is  accelerated  and  at  steady  state  has  the 
Poiseuille  velocity  distribution  given  in  Table  1 

As  shown  in  Table  2,  the  equations  of  continuity  and  motion  are  combined,  and  the  result  is 
transformed  by  a change  of  variables;  the  solution  in  the  form  of  an  infinite  series  is  obtained  by 
separation  of  variables  [3]. 

The  second  example  was  constructed  from  the  first;  the  fluid  was  assumed  to  have  a Poiseuille 
velocity  distribution  initially,  and  the  imposed  pressure  gradient  had  the  same  magnitude  but  opposite 
sign.  The  effect  was  to  decelerate  and  reverse  the  flow.  The  equation  of  motion  and  related  boundary 
conditions  and  initial  conditions  are  listed  in  Table  3. 
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Unsteady— State  Laminar  Velocity  Profiles 
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Unsteady— State  Laminar  Velocity  Profiles 


TABLE  1 

POISEUILLE  VELOCITY  PROFILE 


Dimensional  Form 

Velocity  Profile 

Vz  — Vmax  [1  ~ 

Maximum  Velocity 

r2 

vmax  = (P^  - Pj^)  JJl 

Dimensionless  Form 

Dimensionless  Velocity 

0 — Vz/vmax 

Dimensionless  Position 
Variable 

^ = r/R 

Dimensionless 
Velocity  Profile 

^ = 1 - ^2 
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Unsteady— State  Laminar  Velocity  Profiles 


TABLE  2 

EQUATIONS  OF  CONTINUITY  AND  MOTION 
IN  CYLINDRICAL  COORDINATES  (r,^,z) 
FOR  A FLUID  INITIALLY  AT  REST 


EQUATION  OF  CONTINUITY  (vr,v ^ = 0,  /?  and  //  constant  ) 

= 0 

Qwrr 

EQUATION  OF  MOTION  (vj,  = 0,  /?  and  /z  constant) 


^Vz  _ ^ . 


fi  d 

> 

\ 

7^ 

L 

TRANSFORMATION  OF  VARIABLES 


0 = 


Vz 

Vmax 


0 


e = r/R 


at 


TRANSFORMATION  OF  EQUATION  OF  MOTION 

1 d 


/If. 


or 


1# 


1 d(t)  , 


INITIAL  BOUNDARY  CONDITIONS 


I.C.: 

r = 0 

o 

II 

B.C.  1: 

II 

o 

B.C.  2: 

1 = 1 

o 

II 

ANALYTICAL  SOLUTION  [3] 


00 


n— 1 


where  Qn  are  the  roots  of  J (o;) 

0 
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Unsteady— State  Laminar  Velocity  Profiles 


TABLE  3 


EQUATIONS  OF  CONTINUITY  AND  MOTION 
IN  CYLINDRICAL  COORDINATES  (r,^,z) 
FOR  A FLUID  INITIALLY  IN  MOTION 


EQUATION  OF  CONTINUITY  (vr,v^  = 0,  p and  p constant) 

d 

-g^P^z)  = 0 
gr 

EQUATION  OF  MOTION  (vj-,  = 0»  and  p constant) 


5vz  _ ^ , n 


i d 

r ~Ui 

L J J 

TRANSFORMATION  OF  VARIABLES 


<j>  = 


vz 

Vmax 


e = r/R 


pi 


TRANSFORMATION  OF  EQUATION  OF  MOTION 

1 d 


4 + 


or 


INITIAL  BOUNDARY  CONDITIONS 


I.C.; 

r = 

= 0 

o 

II 

B.C.  1: 

= 0 

o 

II 

B.C.  2: 

= 1 

o 

II 
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Unsteady— State  Laminar  Velocity  Profiles 


COMPARISON  OF  NUMERICAL  AND  ANALYTICAL  SOLUTIONS 


Velocity  profiles,  Figure  2,  for  the  example  in  which  the  fluid  was  initially  at  rest, 
calculated  on  the  interval  0 to  1 with  steps  of  0.1,  [0  (0.1)  1.0],  for  the  equation  of  motion 
form 


to  eliminate  numerical  difficulties  at  ^ = 0.  Results  from  HERCOL  are  denoted  by  asterisks 

analytical  solutions  are  shown  as  dashed  curves  ( ),  and  the  steady— state  solution  is  shown  as  a 

solid  line  ( ).  Tables  4 and  5 list  values  of  (f)  for  ^ on  the  interval  [0,  (0.1),  1.0]  calculated  from 

HERCOL,  and  from  the  analytical  solution;  for  practical  purposes,  the  results  are  identical. 

Figure  3 is  a plot  of  the  velocity  profiles  for  the  second  example;  the  fluid  is  initially  in  motion 
with  a Poiseuille  velocity  distribution  and  a pressure  gradient  with  the  same  magnitude  as  in  the  first 
example  but  with  opposite  sign  as  imposed  on  the  flow  at  T = 0.  At  r = 0.05,  the  velocity  profile  has 
been  distorted  and  at  r = 0.20,  the  flow  is  in  the  opposite  direction.  At  steady— state,  the  velocity 
profile  is  symmetric  about  ^ = 0 with  respect  to  the  initial  velocity  profile.  A condition  such  as  this 
may  occur  in  the  discharge  lines  of  a surge— type  gas  compressor. 


were 

written  in  the 
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Unsteady— State  Laminar  Velocity  Profiles 


Figure  2.  Velocity  profiles  for  unsteady  lonninar  flow 
in  a circular  tube,  fluid  at  rest  initiallyo 
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Unsteady— State  Laminar  Velocity  Profiles 


Figure  3.  Velocity  profiles  for  unsteady  laminar  flow 

in  a cylindrical  tube,  fluid  in  motion  initially. 
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Unsteady — State  Laminar  Velocity  Profiles 


TABLE  4 

COMPARISON  OF  NUMERICAL  AND  ANALYTICAL  SOLUTIONS  OF 


ON  THE  INTERVAL  [0  (0.1)  1.0]  AT  r = 0.05 


^ HERCOL 

^ analytical 

% Difference 

0.00 

0.19962 

0.19959 

0.015 

0.10 

0.19948 

0.19946 

0.010 

0.20 

0.19897 

0.19894 

0.015 

0.30 

0.19772 

0.19769 

0.015 

0.40 

0.19497 

0.19495 

0.010 

0.50 

0.18935 

0.18933 

0.011 

0.60 

0.17861 

0.17859 

0.011 

0.70 

0.15936 

0.15934 

0.013 

0.80 

0.12701 

0.12698 

0.024 

0.90 

0.07594 

0.07592 

0.026 

1.00 

0.00000 

0.00000 

HERCOL 


: numerical  results  from  HERCOL 


<i>  , • analytical  solution 

^ analytical 
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Unsteady— State  Laminar  Velocity  Profiles 


TABLE  5 

COMPARISON  OF  NUMERICAL  AND  ANALYTICAL  SOLUTIONS  OF 


ON  THE  INTERVAL  [0  (0.1)  1.0]  AT  r = 0.20 


e 

^ HERCOL 

^ analytical 

% Difference 

0.00 

0.65181 

0.65178 

0.005 

0.10 

0.64680 

0.64678 

0.003 

0.20 

0.63158 

0.63156 

0.003 

0.30 

0.60551 

0.60549 

0.003 

0.40 

0.56759 

0.56758 

0.002 

0.50 

0.51646 

0.51645 

0.002 

0.60 

0.45049 

0.45048 

0.002 

0.70 

0.36782 

0.36781 

0.003 

0.80 

0.26650 

0.26650 

0 

0.90 

0.14454 

0.14453 

0.007 

1.00 

0.00000 

0.00000 

0 


HERCOL 


: numerical  results  from  HERCOL 


6 , • analytical  solution 

^ analytical 
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Unsteady— State  Laminar  Velocity  Profiles 


CONCLUSIONS 


The  performance  of  HERCOL  was  quite  satisfactory  for  the  solution  of  the  equation  of  motion  in 
cylindrical  coordinates, 


with  combined  Dirichlet,  (f) 


, and  Neumann 


•If 


, boundary  conditions. 


The  code  will  be  useful  for  generation  of  solutions  for  similar  problems  in  which  the  initial  and  boundary 
conditions  preclude  analytical  solutions,  the  equation  of  state  is  not  a simple  function,  or  the  transport 
properties  of  the  fluid  are  not  constant. 
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Unsteady— State  Laminar  Velocity  Profiles 

NOMENCLATURE 

Symbols 


p = 

= pressure 

r = 

= radial  coordinate,  L 

R = 

= tube  radius,  L 

t = 

= time,  t 

V = 

= velocity,  L/t 

z = 

= rectangular  coordinate,  L 

a - 

= roots  of  Bessel  function  of  the  first  kind,  zero  order 

P = 

= density,  M/L^ 

= dimensionless  position  variable 

<!>-- 

= dimensionless  volocity  variable 

= azmithual  coordinate,  radians 

T = 

= dimensionless  time  variable 

( Subscripts 

analytical  analytical  solution 

HERCOL  numerical  solution 


L 

station  "L" 

max 

maximum 

z 

z— component 

0 

station  "0" 

Mathematical  Functions 

Jo  = Bessel  function  of  the  first  kind,  zero  order 
Jl  = Bessel  function  of  the  first  kind,  first  order 
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LISTING  OF  SOURCE  CODE  FOR  CALCULATION  OF 
UNSTEADY  LAMINAR  VELOCITY  PROFILES 
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C Name: 


BSL126.FOR 


C Required: 


HERCOL.FOR 


Purpose: 


Solve  for  unsteady  laminar  flow  in  a 
circular  duct.  Example  4.1-1  BSL 


Keywords : 


fluid  flow,  laminar,  duct,  unsteady 


Type: 


Program 


Status : Experimental 

Reference:  Bird,  R.B.,  W.E.  Stewart,  and  E.N.  Lightfoot, 

"Transport  Phenomena",  John  Wiley  & Sons,  New  York 
(1960) , p.  126 


C Version:  011591 


PROGRAM  HCTEST 


PARAMETER  (NU1=81,  NU2=4,  LW=20) 

REAL  RWORK(IOOOO) , XPTS(NUl),  U(NU1,  NU2) 

REAL  WL(LW),  WR(LW) 

REAL  XO(NUl),  UX(NU1,  NU2) , UERR(NU1,3),  UMX,  SECOND 
REAL  RTOL,  ATOL,  T,  TOUT,  CPUl,  CPU2 , PI2,  ALPHA,  EPSA,  EPSR 
INTEGER  IWORK(IOOO),  INFO(15),  NODE(2),  LUNIT,  NPRT,  NSYS 
INTEGER  NPTS,  LRW,  LIW,  MESS,  L,  NPDE,  M,  K,  IDID,  NODE2 , I 

COMMON  /FCOM/  NSYS 

COMMON  /INI COM/  ALPHA 

COMMON  /TOLCOM/  NODE2 (2 ) , EPSR, EPSA 

C INTEGER  NU1,NU2,LW 

DATA  LRW,  LIW/20000,  1000/ 

C 

OPEN ( 9 , FILE= ’ BSERRM. OUT ' , STATUS= ' UNKNOWN ' ) 

OPEN ( 10 , FILE= ' BSCALCS . OUT ' , STATUS= ' UNKNOWN ' ) 

C DO  600,  MM  = 1,5 

MESS  = 10 


C WRITE  ERROR  MESSAGES  ON  UNIT  9 
LUNIT  = 9 
NSYS  = 2 


C INITIAL  VALUE  OF  THE  TIME-LIKE  VARIABLE 
T = 0. 
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C FINAL  VALUE  OF  THE  INTEGRATION  VARIABLE 


TOUT  =0.20 

C INFO  IS  USED  TO  SET  OPTIONS 

DO  005,  L = 1,15 
INFO(L)  = 0 
005  CONTINUE 

ALPHA  = 10** (MM  + 1) 

WRITE  (MESS,  200)  ALPHA 

C NUMBER  OF  ODE  VARIABLES  WL  AT  LEFT  BOUNDARY  XPTS(l) 

NODE(l)  = 0 

C NUMBER  OF  ODE  VARIABLES  WR  AT  RIGHT  BOUNDARY  XPTS(NPTS) 
NODE (2)  = 0 

C RELATIVE  AND  ABSOLUTE  ERROR  TOLERANCES 

RTOL=1.0E-4 

ATOL=1.0E-4 

C NUMBER  OF  COMPONENTS  IN  THE  VECTOR  OF  UNKNOWNS,  U 
NPDE  = 1 

C NUMBER  OF  BREAKPOINTS,  XPTS ( 1) . . . XPTS (NPTS) 

NPTS  =11 

C SET  UP  THE  MESH  AND  SPECIFY  THE  INITIAL  VALUE  OF  U 
DO  010  1=1, NPTS 

XPTS (I)  = REAL (I-l) /REAL (NPTS-1) 

C Bird,  Stewart,  and  Lightfoot  (4.1-28) 

U(I,1)  =0.0 

C Bird,  Stewart,  and  Lightfoot  (4.1-27) 

UX(I,1)  = 0.0 

010  CONTINUE 

EPSR=RTOL 

EPSA=ATOL 
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N0DE2 (1)=N0DE(1) 

N0DE2 (2)=NODE(2) 

C WRITE  OUT  THE  INPUT  PARAMETERS 
WRITE  (MESS, 198) 

WRITE  (MESS,  210)  NPDE,  NPTS,  T,  TOUT,  RTOL,  ATOL,  NSYS , 

$ (INFO(L),  L=l,  15) 

WRITE  (MESS,  220)  NODE 

C TIME  THE  HERCOL  SUBROUTINE 

CPUl  = SECOND 0 

CALL  HERCOL ( INFO , IDID , U , UX , NUl , NPTS , NPDE , WL , WR , NODE , XPTS , 

1 T , TOUT , RTOL , ATOL , RWORK , LRW , I WORK , LIW , LUNIT ) 

CPU2  = SECOND  0 - CPUl 

IF (IDID  .LT.  0)  THEN 

WRITE  (MESS, 230)  IDID,T 
GOTO  600 
END  IF 

WRITE  (MESS,  240)  T,  IDID,  CPU2 

WRITE  (MESS,  250)  RWORK (7 ) , (IWORK(L) , L=ll,  15),  IWORK(8) 


C CALCULATE  THE  ERROR 
IFLAG  = 0 

UMX  =0.0 
SQERR  =0.0 

C COMPUTE  THE  SQUARE  ERROR 
C DO  020  I = 1,NPTS 

C VAL  = THETA(XPTS(I) , T,  IFLAG) 

C UERR(I,1)  = ABS(U(I,1)  ~ VAL) 

C UMX  = MAX(UMX,ABS(UERR(I,1) ) ) 

C SQERR  = SQERR  + UERR(I,1)**2 

020  CONTINUE 

WRITE  (MESS, 370)  NSYS,T,UMX 

WRITE  (MESS, 398) 

DO  030  I = 1,NPTS 


WRITE  (MESS, 405)  I,  XPTS(I),  U(I,NPDE) 

C WRITE  (MESS, 405)  I,  XPTS(I),  U(I,NPDE),  UX(I,NPDE),  UERR(I,NPDE) 

030  CONTINUE 
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C WRITE  (MESS, 201)  SQERR 

600  CONTINUE 

CLOSE (9) 

CLOSE (10) 

C FORMAT  STATEMENTS 

198  FORMAT  ( IX, '++++++++++++++++++++++++' ) 

200  FORMAT  (IX, 'ALPHA  = ’,1PE12.4) 

201  FORMAT  (IX,'  SQUARE  ERROR  = ',1PE12.4) 

210  FORMAT  (/ ' EX  4.1-2  BSL,  P.  126  '/IX, 'NPDE= ' , 12 , ' NPTS= ' , 

$ 13,'  T=',  1PE10.3,'  TOUT=',  ElO . 3/lX, ' RTOL= ' , E9 . 2 , ' ATOL= ' , 

$ E9.2,  ' NSYS=',  I2/1X, 'INFO=' , 1514) 

220  FORMAT  ('  NODE=',  215) 

230  FORMAT (/'  *****  HERCOL  FAILED  IDID=',I6,'  T=',1PE10.3) 

240  FORMAT  ('  T= ' , IP  ElO. 3,'  IDID=',  15,'  CPU= ' , E9.2) 

250  FORMAT  (3X, ' H= ' , IP  E9 . 2 , ' NSTEP= ' , 14,'  NFE= ' , 14,'  NJE=',  14, 
$ ' NEF=',  13,'  NCF=',  13,'  NQ= ' , II) 

370  FORMAT(5X, 'NSYS=' ,12, ' T= ' , IPEIO . 3 , ' MAX  ABS  ERROR= ' , ElO . 3 ) 

398  FORMAT  (IX, 'TABLE  OF  CALCULATED  RESULTS') 

405  FORMAT  ( IX, 14 , 2X , 5 ( 1PE12 . 4 ) ) 

610  FORMATC  UERR(M,I)  M= ' , 12/ ( 5X,  1P8E9 . 2 ) ) 

620  FORMATC  UERR(I)  '/ (5X,  1P8E9 . 2 ) ) 

END 


C 
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SUBROUTINE  FUN(T,  XC,  UT,  U,  UX,  UXX,  NCPT,  NPDE,  FT,  IRES) 

C FOR  THE  HCTEST  OF  HERCOL 

REAL  T,XC(NCPT),  U(NCPT,  NPDE) , UX(NCPT,  NPDE) , 

$ UXX(NCPT,  NPDE),  FT(NCPT,  NPDE),  UT(NCPT,  NPDE) 

COMMON  /FCOM/  NSYS 

INTEGER  NCPT,  NPDE,  IRES 
INTEGER  M,  I 
INTEGER  NSYS 

DO  020  M = 1,  NPDE 

DO  010  1=1,  NCPT 

C Bird,  Stewart,  and  Light foot  (4.1-21) 

FT(I,  M)  = XC(I) *UT(I,M)  - XC(I)  * UXX(I,M)  - UX(I,M)  - 4 * XC(I) 

010  CONTINUE 

020  CONTINUE 

RETURN 

END 


SUBROUTINE  BDYLFT(T,  UT,  U,  UX,  NPDE,  WT,  W,  B,  IRES) 

FOR  THE  HCTEST  OF  HERCOL 

REAL  T,  UT(NPDE),  U(NPDE),  UX(NPDE),  WT(*),  W(*) , B(*),  PI2 

COMMON  /FCOM/  NSYS 

INTEGER  NPDE,  IRES 
INTEGER  NSYS 

C FINITE  VALUE  OF  U AT  LEFT-HAND  BOUNDARY 

B(l)  = UX(1) 

RETURN 

END 

C ======================================================= 

SUBROUTINE  BDYRHT(T,  UT,  U,  UX,  NPDE,  WT,  W,  B,  IRES) 


20 


REAL  T,  UT(NPDE),  U(NPDE),  UX(NPDE) , WT(*),  W(*) , B(*),  PI2 
INTEGER  NPDE , IRES 
INTEGER  NSYS 

COMMON  /FCOM/  NSYS 
COMMON  /INICOM/  ALPHA 

C U = 0 AT  RIGHT-HAND  BOUNDARY 

B(l)  = U(l) 

RETURN 

END 
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</></> -y> -c/> 


C 


SUBROUTINE  JACOB (T,  XC,  UT,  U,  UX,  UXX,  WTL,  WL,  WTR,  WR, 

$ NCPT,  NPDE,  DFDU,  DBLDW , NBL , NWL , DBLDU,  DBRDW,  NBR,NWR, 

1 DBRDU) 

C PROVIDE  AN  ANALYTIC  JACOBIAN  FOR  THE  HCTEST  CODE 

REAL  XC(NCPT),  UT(NCPT,  NPDE),  U(NCPT,  NPDE), 

UX(NCPT,  NPDE),  UXX(NCPT,  NPDE),  WTL(*) , WL(*) , WTR(*) , 
WR(*),  DFDU(NCPT,  NPDE,  NPDE,  4),  DBLDW(NBL,  NWL,  2), 
DBLDU(NBL,  NPDE,  3),  DBRDW(NBR,  NWR,  2), 

DBRDU (NBR,  NPDE,  3) , T 

COMMON  /FCOM/  NSYS 

INTEGER  NCPT,  NPDE, NBL, NBR, NWL, NWR 


INTEGER 

K,M 

INTEGER 

NSYS 

DBLDU ( Iv 

1, 

1)  = 

l.OEO 

DBLDU ( 1 , 

1/ 

2)  ^ 

-l.OEO 

DBRDU ( 1 , 

1)  = 

l.OEO 

DBRDU (1, 

1, 

2)  = 

-l.OEO 

RETURN 

END 
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SUBROUTINE  SETTOL ( RTOLK , ATOLK , RTOLW , ATOLW , NPTS , NPDE ) 

C SET  ERROR  TOLERANCE  ARRAY  FOR  TEST  CASES 

REAL  RTOLK ( NPTS , NPDE ) , ATOLK ( NPTS , NPDE ) , 

1 RTOLW(*),  ATOLW(*) 

REAL  EPSR,EPSA 

COMMON  /TOLCOM/  NODE ( 2 ) , EPSR, EPSA 
COMMON  /FCOM/  NSYS 

INTEGER  NDl, NPDE, NPTS, K,M 
INTEGER  NSYS, NODE 
C 

DO  20  K=1,NPTS 
DO  20  M=1,NPDE 

RTOLK (K,M)=EPSR 
ATOLK(K,M)=EPSA 

20  CONTINUE 


IF(NODE(l)  .NE.  0)  THEN 

DO  30  K=l,NODE(l) 

RTOLW (K)=EPSR 
30  ATOLW(K)=EPSA 

END  IF 

IF (NODE (2)  .GT.  0)THEN 
ND1=N0DE(1) 

DO  40  K=ND1+1,N0DE(2)+ND1 

RTOLW (K)=EPSR 
ATOLW (K)=EPSA 

40  CONTINUE 

END  IF 
RETURN 

END 
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SUBROUTINE  UINIT ( XC , UC , NCPT , NPDE , WL , WR) 

REAL  XC(NCPT) ,UC (NCPT, NPDE) ,WL(*) ,WR(*) ,PI2 
REAL  ALPHA 

COMMON  /INICOM/  ALPHA 
COMMON  /FCOM/  NSYS 

INTEGER  NCPT, NPDE, K,M 
INTEGER  NSYS 

DO  100  K = 1,  NCPT 

UC(K,  1)  = XC(K) **2 

100  CONTINUE 

RETURN 

END 


C 


END  OF  TEST  ROUTINE 
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